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Abstract. We formulate a theoretical framework to describe multiparticle current transport in planar superconducting 
tunnel junctions with diffusive electrodes. The approach is based on direct solving of quasiclassical Keldysh-Green 
function equations for nonequilibrium superconductors, and consists of a combination of a circuit theory analysis and 
improved perturbation expansion. The theory predicts much greater scaling parameter for the subharmonic gap structure 
of the tunnel current in diffusive junctions compared to the one in ballistic junctions and mesoscopic constrictions with 
the same barrier transparency. 
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Figure 1. One-dimensional (a) and planar (b) models of the tunnel junction. 



1. Introduction 

Multiparticle tunnelling (MPT) is known to be a mechanism of dissipative current transport in superconducting 
tunnel junctions at the subgap applied voltage eV < 2A, and at small temperature T <C A (TJ. It has been shown 
in J2] |3] that the MPT is completely equivalent to the coherent multiple Andreev reflection mechanism (MAR) 
||4] |5| • Here we use the term MPT to emphasize the low transparency, tunnel junction limit, leaving the term 
MAR for a general case of transparent weak links. Each MPT event can be considered as a chain of multiple 
Andreev reflections (6j accompanied by the transfer of n electrons through the tunnel barrier and eventually results 
in the creation of two quasiparticles which contribute to the dissipative current. This process manifests itself in the 
set of the current steps at eV = 2A/n - the subharmonic gap structure (SGS) - which is commonly observed in 
planar junctions with tunnel barriers [7, 8, 9, 10], and in tunable mesoscopic constrictions in the tunnelling regime 
ifTTl [121 n~3l . The theory predicts the scaling D/2 between the neighboring current steps in the tunnelling SGS, 
where D is the bare transparency of the junction tunnel barrier JT] 12 [3] IU . However, in the experiment, the SGS 
scaling parameter is generally much larger. A common explanation for this enhancement is the imperfection of the 
junction insulating layer [6 8, 14] . In our previous paper lfT31l the attention was drawn to an alternative explanation: 
effect of disorder in planar junction electrodes; it was predicted that the SGS scaling parameter for planar diffusive 
junctions is enhanced by a factor ~ ^o/£, or even ~ E,q /id depending on the junction geometry (^o is the coherence 
length, i is the elastic mean free path, d is the thickness of electrodes). 

In this paper we present a detailed theory of the MPT in planar Josephson tunnel junctions with diffusive 
thin-film electrodes and extensively discuss the details of the behaviour of the MPT currents and the relevant 
asymptotical methods. The theory is based on the direct solving of the diffusive equations of nonequilibrium 
superconductivity ifTBI . The main difficulty with this approach is related to essentially nonstationary character of 
the Josephson tunnel transport. While analytical and numerical methods are well developed for solving stationary 
Usadel equation [17| and nonequilibrium Keldysh-Usadel equation [18], the nonstationary problem was so far 
studied only numerically [19]. The central model assumption made in this paper and relevant for the tunnel 
regime concerns discrimination of nonzero harmonics of the Green's functions and the distribution function. This 
approximation turns the originally difficult problem into an analytically tractable one, and at the same time it 
captures all qualitative, and to large extent quantitative properties of the tunnelling SGS. Within this approximation 
we are able to develop a relatively simple and physically appealing calculation scheme, which combines the 
improved iterative procedure for evaluating the tunnelling density of states (DOS), and the circuit theory methods 
lfT8ll20 j for evaluating the dc current. 

The resulting physical picture of the MPT is as follows: The tunnelling processes create the local nonzero 
DOS inside the bulk energy gap in the vicinity of the tunnel junction. This allows quasiparticles to overcome 
the energy gap at a small applied voltage in several steps, by repeated bouncing between the junction electrodes 
(MAR). The spectral current through the energy gap, which determines the net charge current, is calculated by 
considering an effective circuit theory network representing the tunnelling process. 

The effect of substantial enhancement of the SGS scaling factor can be qualitatively understood from 
a mesoscopic picture of the diffusive tunnelling transport, namely, tunnelling through a set of independent 
quasi-ballistic conducting channels with randomly distributed transparencies l20l l2"fl . Within this picture, the 
contribution of each channel can be evaluated using the ballistic MAR theory J21[3j|4l. In constrictions with 
length L ^> £ the transparencies are spread over the interval ~ (L/£)D ^> D ll22l . which implies that the junction 
transparency, and hence the SGS scaling factor, are effectively enhanced by the factor Ljl. This explanation, 
however, is valid only for short constrictions, L <C ^o, while it does not apply to planar tunnel junctions with 
overlapping thin-film electrodes commonly used in experiments and shown in figureQ] In these structures, massive 



Multiparticle tunnelling in diffusive superconducting junctions 



3 



pads (reservoirs) are fairly far from the junction (L ^> in such a situation the effective junction length is 
defined by the scale of spatial variation of the Green's function, i.e., by the coherence length. If the size Ly of the 
overlapping parts of the electrodes is comparable to the electrode thickness d, then the junction can be considered 
as an effectively one-dimensional (ID) one (see figure [Tfa)); in this case, the SGS enhancement factor becomes 
£,0 /£. If the junction cross-section is much larger than the cross-section of the electrodes (see figure |TJb)), the 
current concentrates not at the junction but rather in electrodes |23|, and an additional enhancement factor %o/d 
appears in the SGS scaling, which coincides with the result of rigorous calculation in this paper. Remarkably, 
this enhancement concerns only the multiparticle currents, while the single-particle current is not affected and is 
proportional to the bare transparency D l24l . 

The paper is organized as follows. In its major part we develop a theory for the ID junctions, figureQJa), and 
discuss the extension to the planar junctions, figure 02b), towards the end, in section |U We start with discussion 
of basic equations and adopted approximations in section|2] Then we construct the circuit theory in section[3j and 
develop a perturbation theory for the DOS in section|5] The MPT currents are calculated in sections [6] and [7J the 
latter section includes also the calculation of the excess current. The effect of neglected harmonics in the Keldysh 
and Green's functions is evaluated in section[8] In section|9]we discuss the results and possible implications of the 
theory. 



2. ID junction model 



2.1. Basic equations 

The model of tunnel junction we are first going to study is depicted in figure [Ha) and consists of a tunnel barrier 
with the transparency D attached to bulk superconducting electrodes via two superconducting leads (— L < x < 
and < x < L). We will consider a diffusive limit, in which the elastic scattering length t is much smaller than 
the coherence length = \/33/2A, where T> is the diffusion coefficient (we assume h = Icq = 1). We assume the 
length L of the leads to be much larger than <^o, and their width to be much smaller than the Josephson penetration 
depth which implies homogeneity of the current along the junction. Similar model has been considered in [25 1 in 
study of the dc Josephson effect in tunnel structures. 

Under these conditions, the microscopic calculation of the electric current 7(f) requires solution of the 
diffusive equations of nonequilibrium superconductivity |16| for the 4x4 matrix two-time Keldysh-Green's 
function G(x,t\ ,t 2 ) in the leads, 

[H,oG\=YDdJ, J = God x G, GoG = S(h-t 2 ), (1) 

H= [ia z d ti -e(p+A(h)]8(h-t 2 ), A = j a *Ho y A, 

where <p is the electric potential, A and (j) are the modulus and the phase of the order parameter, respectively, <7,- are 
the Pauli matrices, d x denotes partial derivative over the variable x and 

3 R G K 







G K =g R of-fog A . (2) 



Here g R ' A are the 2x2 Nambu matrix retarded and advanced Green's functions, and f — f + <7,/_ is the matrix 
distribution function (we use 'check' for 4x4 and 'hat' for 2x2 matrices). The multiplication procedure in (fl} 
and (O involves the time convolution 

(AoB)(t u t 2 )= [ A(t u t)B(t,t 2 )dt. (3) 

For arbitrary tunnel barrier, the function G and the matrix current J at the left (x = —0) and the right (x = +0) 
sides of the tunnel junction are connected via the generalized boundary condition by Nazarov [20|, 

_ _J_ /-i Dp(D)dD[G_ 0; oG +0 ] 
J -° 7+0 2g N Rj 1 + £({G_o,oG +0 }-2)' 1 

where p (D) is the distribution of the transparencies of the conducting channels of the barrier (/J Dp (D) dD =1), 
R is the junction resistance and gN is the normal conductivity of the leads per unit length. Assuming the absence of 
high-transparent channels with D ~ 1 and considering p (D) to be localized around the small value of D <C 1 (tunnel 
limit), we can neglect the anti-commutator term in (@), thus arriving to the Kupriyanov-Lukichev's boundary 
condition [26|, 

/-o-/+o=(2^)" 1 [G_o,oG +0 ]. (5) 
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The electric current is related to the Keldysh component of the matrix current /as I(t) — (ng^/4-e)Tra z J K (x,t,t), 
and thus it can be expressed through the boundary value Jq, 

/(f) = (7r/8 e /?)Trcj,[G_ ,oG +0 ] K (f,f). (6) 
Equations (Q]) can be decomposed into the diffusion equations for the Green's functions, 

[H, og] = iHdJ, J = god x g, g o | = 8 (f i - 1 2 ) , (7) 
and the equation for the Keldysh component G K , 

[H,oG K ] ^iVdJ*, J K =g R od x G K + G K od x g A , (8a) 

| R oG K + G K o| A = 0. (8b) 
The boundary conditions for the functions g and G K at the tunnel barrier follow from ©, 

f = (W/$o)[g-o,°g+o\, (9a) 

^-(W/^o)[G-o,oG +0 ] K . (9b) 
In (|9oT > and d9fct . the transparency parameter W is defined as 

W=.R(§,)/2K=(3§>/4*)I>»Z>, (10) 

where R(%o) = 4o/gN is the normal resistance of the piece of the lead with the length <^o- It has been shown in ll25l 
that it is the parameter W rather than the barrier transparency D that plays the role of a true transparency parameter 
in diffusive tunnel junctions. We will consider the limit W <C 1, which corresponds to the conventional tunnelling 
concept. In this case, according to d9ob and the gradients of all functions are small. Within the tunnel model, 
which assumes W to be the smallest parameter in the theory, these gradients are neglected, and the functions g and 
/ are taken local-equilibrium within the leads. In our consideration, we will lift this assumption and suppose the 
local-equilibrium form of these functions only within the bulk electrodes (reservoirs). Attributing the reference 
point for the phase, (j) — 0, to the left electrode, x = —L, these functions in the right electrode, x = L, are given by 
relations 

g{E,t) = a z u(E + a z eV)+ie ia ^ t) a y v(E), (11a) 
{u,v) = (E,A)/^, § R > A = p±iO) 2 -A 2 ] 1/2 , (lib) 
f(E)=Umh[(E + o z eV)/2T], (11c) 
written in terms of the mixed Wigner representation A (E,t) of the two-time functions, 

AE 



where the variable E has the meaning of the quasiparticle energy, and t = (t\ +?2)/2 is a real time. Similar 
equations, with 0=0 and V = 0, apply to the left electrode, x = —L. 

Because of the small value of the tunnelling parameter W one can neglect variations of the superconducting 
phase along the leads, as well as the charge imbalance function /_ proportional to a small electric field penetrating 
the superconducting leads. Furthermore, the small value of the superfluid momentum in the superconducting 
leads, p s ~ W ll25l . allows us to neglect a small effect of the energy gap suppression by the superfluid momentum 
(~ p\^ ~ W 4 / 3 ll27l ). Within such an approximation, the coefficients in ([U at the left lead, x < 0, are time- 
independent functions. At x > 0, applying the gauge transformation [28 1 G(fi,?2) = 5^(?i)G(fi,f2)S(f2) with a 
unitary operator 5(f) = exp[i<7 z 0(f)/2] to the function G, we exclude the time-dependent phase and the electric 

potential from the equations for the function G and corresponding boundary conditions at x = L, which then 
become similar to the equations for G(x) at x < and the boundary conditions at x = —L. This results in 

the symmetry relation G(x) = G(—x), which allows us to replace the function G+o in the boundary condition 
<(5j and in the expression © for the electric current by the inverse gauge transformation of the function G_o, 

G+o (Lo = 5(f 1 )6 +0 5 t (f 2 ) = 5(f 1 )G_o5 t (f 2 ). 

As the result, the problem is reduced to the solution of a static equation within the left lead for the function 
G(x,fi,f2), completed with the time-dependent boundary condition (O at the tunnel barrier. Similar approach is 
used in the theory of ballistic point contacts [2 | where the Josephson coupling is described by an effective time- 
dependent matching condition for the gauge-transformed Bogolyubov-de Gennes equations in the leads. 
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It is convenient to expand all functions over harmonics of the Josephson frequency, A(E,t) = 
Y im A(E,m)exp(—2ieVmt), using the following rules for representation of the products and gauge-transformed 
values, 

(AoB)(E,m) (12) 
= £ m , A [E + eV (m - m') , m'] B{E- eVm' , m - m') , 

\ A 2 i(E,m-l), A 2 2{E-eV,m) J 
In such a representation, the expression © for the dc current / has the following form 

AJ7 — 

^ Tr(/!oG K -/zoG K )(£,0) 



_» 16eR 

^Y^ m [h{E,m)W{E,-m) ( 14 > 

-Ti(E,m)G K (E,-m)] , h = <j z f - g A a z , 

where all functions are taken at the boundary x = — 0. We will adopt the same convention in the most of following 
equations and assume the spatial coordinate to be taken at the boundary. 

2.2. Zero-harmonic model 

Solving a system of nonlinear differential equations (|7Ti-(|%ll generally can be fulfilled only numerically even in 
the ID case. The analytical solutions can be constructed in the adiabatic limit of small applied voltage eV <C A 
||29l . To make the problem tractable at larger voltages eV ~ A, we make use of the observation that the amplitudes 
of high-order harmonics of the function G are small in the tunnelling limit W -C 1: the amplitude of the mth 
harmonic decreases with m as W m . This suggests that zero harmonics m = Q play the key role in (TBI) , while the 
high-order harmonics are neglected. Thus we adopt an approximation scheme, in which only the zero harmonics 
of the functions g and G K are kept. It turns out that such an approximation is sufficiently powerful to recover all 
specific features of the MPT currents, and to give satisfactory description of the SGS of the net tunnel current. 
Furthermore, our analysis of the correction due to the first harmonics in section [8] shows that the zero-harmonic 
model may give rather good quantitative agreement with the result of full numerical calculation. 
For the matrix structure of the zero harmonic of the function g in ( TBI ), we adopt the form 

g(E,0)=o z u(E)+io y v(E), u 2 -v 2 = l, (15) 

which is similar to the Green's function structure ( II lab in the left electrode, though u and v differ from their 
equilibrium values in ( 11 1M . It is possible to prove, using the normalization condition in (O, that the zero harmonic 
of matrix g is traceless, and its a r -component is much smaller (at least by W 2 ) than the "main" components u and 
V, which have zero order in the parameter W. Within the same approximation, the Keldysh function has the form 
G K (£,0) = 2f{a z N + ia y M) , where N(E) = Rew R is the density of states (DOS) normalized to its value in the 
normal state and M(E) = Rev R . In what follows, we will express the advanced functions through the retarded 
ones, (h,v) a = — (m,v) r *, using the relation g A — —G z g R1 <G z , and omit the superscript R, assuming all Green's 
functions to be retarded. 

Retaining only zero harmonics of the functions g and G K in ( TBI ), we find that only the diagonal parts h(E) 
and G K (E) of the corresponding matrices enter the dc current 

dF 

TrY j (l+kcj z )[h(E)G K (E+keV) (16) 
-h(E + keV)G K {E)]. 

By introducing the distribution function n = ^(1 — /) which approaches the Fermi function np in equilibrium, 
equation ( fTol l exactly transforms to the standard form for the tunnel current, 

/= / — N(E)N(E-eV)[n(E-eV)-n(E)], (17) 
J-™ eR 

with that essential difference that the DOS and the distribution function are not given, but are to be computed from 
the Keldysh-Green's function equations. To zero order in the tunnelling parameter, the DOS has the BCS form 
Ns(E) = Re (E /£,), and the distribution function is the equilibrium one, n = np. In this approximation, equation 
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( fT7| i recovers the single-particle current of the tunnel model (24]. At zero temperature this current acquires the 
form 

reV-A 

1= \ —N s {E)N s (E-eV) (18) 
J A eR 

and turns to zero at eV < 2A, having a sharp onset at eV = 2A, I\ (2A) = nA/2eR. 

To calculate the current at smaller, subgap voltages eV < 2A, one has to calculate the tunnelling corrections 
to the BCS DOS, and to find the nonequilibrium distribution function. 

3. Circuit representation of boundary condition 

We start with evaluation of the distribution function and develop a circuit theory approach to derive a general 
analytical equation for the current (fTTI i assuming the DOS to be modified in a close vicinity of the tunnel barrier. 

3.1. Kinetic equation and boundary condition 

Using zero harmonic of (I5ab and i9b\ , we obtain the diffusive kinetic equation and the boundary condition, 

d x (D + d x n)=0, (19) 
D+d x n\ x=Q = (2W /^)^ k=±l NN k {n k -n)\ x=Q1 (20) 

where D + = |(1 + \u\ 2 — |v| 2 ) is the dimensionless diffusion coefficient, and the subscript k denotes the energy 
shift: n k (E) = n(E k ) = n(E +keV). It follows from ( fT9] l that D + d x n = const; this constant value can be found 
from (f20b . At x ^> |o the coefficient D + approaches the BCS form D$ = ®(\E\ — A) (©(*) is the Heaviside 
step function) and exactly turns to D§ at the reservoir, x = —L. This implies that at subgap energies \E\ < A, the 
quasiparticle probability current D + d x n turns to zero along the whole lead, and the distribution function is spatially 
homogeneous. Physically, this manifests complete Andreev reflection in terms of quasiparticle flows in diffusive 
structures. 

Outside the gap, equation ( fT9l has no bound solutions: the distribution function grows linearly with x far 
from the junction. Such a growth is limited in practice by inelastic collisions, which provide relaxation of n(x,E) 
to the equilibrium value np at the distance of inelastic scattering length l £ . To simplify the problem, we consider, 
instead of including a complicated inelastic collision term, a junction geometry with short enough leads having the 
length L^il £ (but still L <jjo) and connected at x = ±L to equilibrium reservoirs. Within this model, the Green's 
functions, which change at x ~ are not affected by the finite length of the leads, while the reservoirs impose 
the equilibrium boundary conditions for the distribution function, «(±L) = np. At the same time we can neglect 
inelastic collisions inside the leads. Obviously, this model describes a qualitative pattern of inelastic relaxation in 
very long channel (L 3> £ e ) with substitution L — > £ £ in our results. 

Generally, spatial variation of n(x) at |£| > A has two scales: linear x-dependence at x ~ L and fast but 
small variations near the junction due to spatial dependence of the Green's function. Neglecting these variations 
within the main approximation in the parameters W and %o/L, we arrive at the relation D + d x n = [«(0) — np]L~ l . 
Substituting it into the boundary condition ( |20l i and accounting for D + d x n = at \E\ < A, we obtain the equation 

&(\E\- A) [n - np) = r^ k=±l NN k (n k -n). (21) 

In this and following equations, the functions are taken at the boundary x = —0. Equation (fJTJi represents a 
recurrence relation between the values of n(E) at the energies shifted by eV; the nonequilibrium parameter r is 
defined as 

r^2LW/E,o=R t ,/R 1 (22) 

where R^ = L/g^ is the normal resistance of one lead. In practice, the tunnel resistance much exceeds R^, and the 
nonequilibrium parameter is small, r <C 1, which implies that at the energies \E\ > A, the distribution function is 
always close to the Fermi function. 

3.2. Circuit theory 

We split the integral in (fTTJ into pieces of length eV, 

'-L%^-C%^- <23) 



J = TI=-oo h , Jk = («*- 1 - n k ) p k 1 , p k 1 = N k N k - 1 . 
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Figure 2. Circuit theory network with the period eV representing charge transport in diffusive tunnel junctions, 
eV = 2.5A. 



In these notations the recurrence relation ( f2TT > reads 

&(\E k \ - A) [n k - n F (E k )} = r(j k - j k+l ). (24) 

A convenient interpretation of equations (l23l l and d24T i in terms of the circuit theory [18] is given by an infinite 
network in the energy space with the period eV , graphically presented in figure [2] The electric current spectral 
density J(E) consists of partial currents j k which flow through the chain of tunnel "resistors" p k connected to 
adjacent nodes of the network having "potentials" n k and n k -\. At \E\ > A, the nodes are also attached to the 
distributed "equilibrium source" n F (E) through equal resistors r. In this representation the recurrence relation (1241) 
has the meaning of "Kirchhoff rules" for partial currents. 

Below we assume the equilibrium quasiparticle distribution n(E) = n F (E) at \E\ > A, neglecting effect of 
small resistors r. In this limit, the currents j k outside the energy gap, \E k \ > A, vanish at zero temperature since 
«F is piecewise constant. At T ^ 0, these currents describe the effect of thermal excitations. The subgap spectral 
current is conserved, j k = /'a = const, k = 1 — N- , . . . ,N+ (as a consequence of (f24l i). and can be easily computed, 

./A = [«F {E-N- ) - "F (E N+ )] p A 1 , 

N ± (E) = Int [(ATE)/eV] + 1, p A (E) = Y^i-n. Pk, 

where the integers ±A^± are the numbers of the nodes outside the gap nearest to the gap edges, Int(x) denotes 
integer part of x, and the quantity Pa(E) has the meaning of the net subgap resistance. The subgap distribution 
function reads 

n(E) = n F (E N+ ) + [n F (E_ N _ ) - n F (E N+ )] PkP A ^ (25) 

The resulting electric current can be now written in a general form, 

eV dE , , r dE np(E) — n F (Ei ) 

(N_+N + )j A + 2 — n > PV X > , (26) 
o eR Ja eR p\ 

valid for arbitrary voltages and temperatures. Here the first term describes the subgap current, and the second term 
- the current of thermal excitations. 

The magnitude of the subgap current is fully determined by the net subgap resistance; the current is blocked 
when this resistance is infinite, p A = °°, which happens when DOS turns to zero. According to d26i >. the amount 
+ N + of the resistors contributing to the subgap resistance gives the amount of electric charge (in units of e) 
transferred during the tunnelling event. Thus, the circuit with one subgap resistor represents the single-particle 
tunnelling, which can exist only at eV > 2A; in this case, equation ( 126*1 ) reduces to ( fT~8b . At eV < 2A, the subgap 
circuit should consist of at least two resistors (two-particle tunnelling). However, for the BCS DOS this current is 
blocked, and to evaluate the current one has to calculate the tunnelling correction to the DOS within the gap by 
solving the Green's function equations. 



4. Planar junctions 

In this section we will discuss the extension of our approach to a more practical case of planar tunnel junction 
sketched in figure [Tib). This 2D case is more complex; however, it is possible to reduce this problem to the ID 
case by formulating effective boundary conditions at the junction following the method suggested by Volkov |23l . 

Let us suppose that the size of the junction Lj exceeds the coherence length, Lj 3> and, simultaneously, 
does not exceed the Josephson penetration depth. Then the function G in the left-hand side (lhs) of the kinetic 
equation [H,oG] = iDVJ is approximately constant within the junction banks (parts of the junction leads of lengths 
Lj beneath and above the insulator). Then, integrating this equation over the volume of the bottom bank, using the 
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boundary conditions (O and denoting the cross-section area of the lead as Sg, the lead thickness as d, and the area 
of the junction as St, we obtain 

S T d[H,oG] =iV{S T (W/$o)[G-o,oG +0 ]-S e Jt}, or 

[H,oG] = 2iA{W[G_ ,oG +0 ] - {^/L T )J e }, (27) 

where ±0 denotes top and bottom side of the barrier, Ji> is the value of the matrix current at the lead cross-section 
adjoining the junction, and the tunnelling parameter W is defined as 

W=W(^ /d) = (3^/4£d)D. (28) 

As soon as |o *C Lj and ^qJ( ~ W, the last term in d27l > can be assumed to be the smallest one and thus neglected. 
However, this is only true for the Green's component of ( 127] ! [23], 

[H,og]=2iAW[g- ,o§ +0 ], (29) 

whereas for the Keldysh component, the diagonal part of the lhs of d27l > turns to zero (we consider only zero 
harmonics), and therefore the boundary condition for the diagonal part of Jf, which is proportional to D + d x f, has 
the form 

jf = (W t /^)[G-Q,oG +0 } K , Wf = W(L T /d)-»W. (30) 

Equation (l30l > is the boundary condition for the distribution function, which is to be used as described in previous 
section. Solving the kinetic equation within the lead and assuming L ^> Lj, we arrive at the equation similar to 
(1241 with the same parameter of nonequilibrium, 

2W f L = 2L R{^) L T = p L T = Lp = R N 
r %o 2% R d S T R d S(R R ' 

where p is the specific conductivity of the leads. 



5. Perturbation theory for the Green's functions 

To calculate the DOS within the next approximation with respect to the parameter W, we solve the Usadel 
equation for the Green's function g following from ||7). Introducing usual parametrization g = ovexp((7 v 0) and the 
dimensionless coordinate z, we arrive at the equation for the spectral angle 9, 

sinh[0(z)-0s]=i0"(z)sinh0 s , z=x/% (31) 

(the prime sign denotes the derivative over z). With exponential accuracy, the solution of (l3TT l at z < can be 
approximated by the formula for a semi-infinite superconducting wire [ 30 1 

tanh[(0(z)-0 s )/4]=tanh[(0(-O)-0s)/4]exp(fe), (32) 

where k~ 1 (E) — Visinh 0s- Equation (l32t describes the decay of perturbations of the spectral functions at distances 
> <^o from the barrier, where the spectral angle approaches its bulk value 0s = arctanh(A/£'). The boundary 
condition for the spectral angle follows from ( |9a| l, 

O' + W sinh (cosh 0i + cosh 0_j) \ z= _ Q = 0. (33) 
Then the boundary value of can be found from the finite-difference equation following from (l33l and ( 1321 ). 

2£sinh[(0 s - 0)/2] = W sinh (cosh B\ + cosh 0_ i ) . (34) 
A similar result for the spectral angle in planar junction banks follows from ( 1291 , 

£ 2 sinh(0 s -0) =W sinh 0(cosh0!+ cosh 9-i). (35) 

In what follows, we will simultaneously discuss both of the junction geometries, using common notation W 
for both transparency parameters W and W and assuming this quantity to be defined by ( TTOb or ( f28l depending on 
context. 
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Figure 3. Numerically computed DOS and typical subgap circuits describing the two-particle current at eV = 1.2A (a), 
the three-particle current at eV = 0.7A (b), and the four-particle current at eV = 0.55A (c), for the tunnelling parameter 

w = io- 3 . 



5.7. Simple perturbation theory 

Due to the presence of the small parameter W in the right-hand side (rhs) of 1341 . one can suggest the following 
perturbation correction for 9 to first order in W, 

6 = flSs- Wit -1 sinh^s (cosh +cosh0 s ,-i), (36) 

and similar for ([35}. This results in the following expression for the DOS within the BCS gap, 

N(E) = Re(coshe)=w|A/£| a (Af s ,i+/V s ,-i), (37) 



|S| = VA2-£2, 




ID junction; 
planar junction. 



Such approximation will be referred to as the result of a simple perturbation theory (SPT). 

As follows from d37b . the tunnelling coupling extends the DOS inside the gap over the distance eV from the 
gap edges, and scales it down by the factor W, as shown in figure [3] It is clear from l37b that with decreasing 
voltage, at eV < A, the gap will open in the DOS (see figure Hfb)), and further iteration of the finite-difference 
equations 134b or (f35T > for the spectral angle is required. As the result of this iterative scheme, the DOS at small 
enough eV acquires a staircase structure in the energy space: two ladders descending from the bulk gap edges inside 
the subgap region; the height of nth step is W", and the width is eV. In the middle of the gap, \E\ < A — (n — l)eV 
(assuming eV < 2A/n), there is a plateau with the height w 2W" (see figures[3la) and[3jc)). While eV decreases, 
the plateau expands until its size becomes equal to 2eV, then a new pair of steps emerges, which happens when n 
is even. We recall that this deformation of the DOS occurs only locally, at the distance x ~ from the junction. 



5.2. Improved perturbation theory 

The subgap DOS in (l37T i possesses singularities at the gap edges \E\ = A, which causes a divergence of the subgap 
current at relatively large voltage, as we will see later. To eliminate the divergence, we need to apply an improved 
perturbation theory (IPT) to (l34l and (l35l l. in which nonlinearity of the recurrence relations is fully taken into 
account. First, we consider an approximation to 1341 . in which the non-singular terms coshQ-H are replaced with 
their BCS values, but we do not suppose that the difference 6 — 9s is small, 

sinh[(0 s -0)/2] =Wk~ 1 g(E,eV)smhd, (38) 
g(E,eV) = ^(coshOs.! + cosh 8 S -i). 



Multiparticle tunnelling in diffusive superconducting junctions 10 

In the vicinity of the point E = A the function g(E,eV) is regular and thus can be approximated by g(A,eV) 
(it is sufficient to consider only positive values of E due to the symmetry of the spectral functions). Within the 
region \E — A| <C A, the spectral angles and 0s are large, therefore we hold only large exponents exp and exp 0s 
in the hyperbolic functions in the rhs of d38l l. and use the asymptotic expression exp 0s ~ [2 A/ (E — A)] 1 / 2 . Then, 
introducing a dimensionless energy variable e and a normalized spectral function z(e), 

e = (A-E)/2A P \ p(eV) = [^VW] I/3 , 
z(e)=ipexp0, exp0 s = (ipVe) -1 



we reduce d38t to a numerical algebraic equation 



z 3 + {zVe-iy = 0. (40) 

The relevant solution z(e) of ( f40b is determined by the requirement for the asymptotic behaviour at e 3> 1 to 
coincide with the energy dependence z(e) sa e~ 1//2 + ie~ 5//4 given by the direct perturbative expansion (f36b . For 
planar geometry, we directly obtain the function z(e) and the scaling parameter p from ( l35l >. 

z(e) = (e-i)- 1 / 2 , P W = »)] 1/2 - (41) 
The resulting DOS in the region \E — A| <C A, 

N{E,eV)^[2p(eV)}- 1 lmz(e), (42) 
approaches a finite value iV(A,eV) ~ W~ b (eV — 2A)~ b l 2 at £ = A, where 

{2/3, ID junction; 
, (43) 
1/2, planar junction. 

However, in the vicinity of the specific voltage value eV — 2A the DOS diverges, and the calculation procedure 
must be further improved. The problem is caused by the fact that at eV = 2A both the energies E and E — eV in 
(l34l or (l35l l approach the gap edges A and —A, respectively. Therefore one must solve the equation not only for 
N(E), but for N(E — eV) as well. To this end we consider the recurrence (l34l for these two energies and replace 
the nonsingular terms cosh 9\ and cosh 0_2 by their BCS values, 

s -0 Wsinh0, 
sinh — - — = — — — (cosh0_i +cosh0s i i), 

q q sinh (44) 

sinh — '■ = (cosh 0s _9 +cosh 0). 

Using again the fact that the spectral angles are large by modulus in the region \E — A\ <C A, we hold only large 
exponents exp0, exp0 s » [2A/(E - A)] 1 / 2 and exp(-0_ 1 ), exp(-0 s _ 1 ) « [2A/(eV - A - E)] l l 2 in the rhs of 
(l44l . Then, introducing dimensionless energy and voltage variables e and £2, and normalized spectral functions 

z(e) andz(e), 

e = (E-A)/2Aq 2 , Q.= (eV-2A)/2Aq 2 , 
exp9=zq~\ exp(-0_i)=z ? - 1 ) q = W 2/5 /2. 

we obtain algebraic equations for the functions z(e) and z(e), which reduce to a single equation for the function 

z(e), 

[z 3 (l-zR)- 2 + iR 2 ] 2 + iz[(l-zR) 2 -4zR\ =0, (46) 



, 2 



where R(e) — yfs + i0 and R(e) = R(£l — e). The function z can be then found as z = (1 — zR) / vi?- For planar 
geometry, we obtain the equations 

(l-z 2 /? 2 ) 2 (z + i^ 2 )+iz 4 = 0, z=(l-z 2 R 2 )/iz 2 , (47) 

and the scaling parameter q = (W/4) 1 / 3 . According to the definition of z(e) in (|43T >. the solutions of (|46*] | and ( f4Tb 
are related to the boundary values of the DOS as 

N(E) = (2q)- 1 Rez, N^i (E) = (2 q y 1 Re g. (48) 

The results of computation of the DOS based on the numerical solution of the recurrences ( f34b and < f35l > and shown 
in figure[3] quantitatively confirm the results of our asymptotic analysis. 



Multiparticle tunnelling in diffusive superconducting junctions 



11 



6. Multiparticle currents: thresholds 

6.1. Two-particle current 

The existence of the subgap states enables quasiparticles to overcome the energy gap at eV < 2A via two steps 
involving intermediate Andreev reflection at the energy \E\ < A. The population n(E) of the intermediate state is 
generally non-equilibrium, because the subgap quasiparticles cannot access the equilibrium electrodes. In terms 
of the circuit approach, the node k = is disconnected from the equilibrium source, and the subgap current flows 
through two resistances po and p\ (two-particle current), see figure [3|a). The corresponding partial currents are 
equal, jo = j\ = [np(E\) — np(£'-i)]/(po + Pi), and their contribution to I(V) is confined to the energy region 
< E < eV — A (a similar contribution at A < E < eV comes from jo and j_ i ), which leads to the following 
expression for the two-particle current, 

4 reV-A Ap 

h = — , A<eV<2A. 

eR Jo po + Pi ' 

Within the SPT approximation for the subgap DOS function N, this is equivalent to equation 

A 
I 



\yj reV-A 

h = — dE 



eR Jo 



£,£■_, I 

(49) 



where £,(E)= [(E + iO) 2 — A 2 ] 1 / 2 according to (II Ibh The current h increases with voltage and diverges at eV = 2A 
which is the result of the mentioned DOS singularity at E = A. 

At eV — A, the two-particle current possesses a threshold. In the vicinity of the threshold, eV = A + Q., 
< Q. <C A, equation ( |49l simplifies, giving the current threshold value 

h{ A)^(\ AE =^2W/ l( 2A). (50) 
V ' eR Jo ^&-E 2 eR V 1 

6.2. Three-particle current 

At eV < A, a minigap opens in the DOS around the zero energy (see figure Ob)), however, since the number of 
subgap resistors increases up to three: p_i, po and p\ (three-particle current), the current across the minigap 
will persist as long as the network period exceeds the minigap size, eV > 2(A — eV), i.e., at eV > 2A/3. 
The corresponding partial currents are equal, j_i = jo = ji = [n-p{E{) — np(£-2)]/(p-i +Po + PiX an d tne u- 
contribution to I(V) is confined to the energy region A — eV < E < 2eV — A, which leads to the following expression 
for the three-particle current at zero temperature, 

3 f 2eV-A dE 

h = — , 2A/3<eV<A. (51) 

eRJA-eV P-i+Po + Pi 

Taking the subgap DOS functions jV and AL i in the SPT approximation (f37l >. we see that the central resistance 
Po ~ W~ 2 is the largest. Retaining only this resistance, we get 

\EiE- 2 \ 



W 2 r2eV-A 

h = — dE 



eR 



A-eV 



(52) 



While approaching the voltage value eV = A, the current It, infinitely increases due to decreasing distance between 
the upper integration limit and the singular point £ = A. Calculating ( 1521 in the vicinity of the threshold eV = 2A/3, 
we obtain 

3ttAW 2 /9\ a 

h (2A/ 3 ) = — J— ^— ( - ) «2W/ 2 (A). (53) 



2eR 



6.3. Four-particle current 



At eV < 2A/3 the network period becomes smaller than the minigap, and the situation resembles the one 
encountered when the voltage decreased below 2A: we have to calculate the next correction to the DOS ~ W 2 . In 
the equation for the four-particle current, 

8 f 2eV - A dE A 2A 

h = — , - < eV < — (54) 

eRJo P-1+P0 + P1+P2 2~3 

(see figure|3c)), the largest resistances are po ~ Pi ~ W -3 , thus p_i ~ p_2 ~ W -1 can be neglected. The functions 
N±i can be obtained from the SPT equations (|37| | at E = E±\, in which small in the rhs has to be neglected: 
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eV/A 



2 3 4 



Figure 4. Multiparticle currents in planar junction numerically computed for the tunnelling parameter W = 10 



N±i = W\A/%±\ \ a Ns±2- To evaluate the function N, we must perform next iteration step, by substituting these 
values of N±i into the SPT equation (l37l > for N. As the result, we obtain 

^ ' |£ 2 £- 2 | 



gw3 fieV-A 



eR Jo 



While approaching the voltage value eV = 2A/3, the current I4 infinitely increases. In the vicinity of the threshold 
eV = A/2 we have 

/ 4 (A/2) = ^-UJ «2W/ 3 (2A/3). (55) 

From these considerations we conclude that the evaluation of 2n- and (2n + 1 ) -particle currents requires DOS 
recurrences of nth order. As long as the applied voltage decreases below A/n, a new minigap opens in the DOS 
(see discussion in section IBTTT i. and the recurrent procedure should be repeated again. 



7. Multiparticle currents: large voltage 

It follows from the previous section that the multiparticle currents calculated within the SPT approach have finite 
values in the vicinity of their thresholds, but they diverge at the next gap subharmonics: the two-particle current 
diverges at eV = 2A, the three-particle current diverges at eV = A, and so on. These divergences are caused by 
singularity of the SPT correction to the tunnelling DOS at the point E = A which enters the integration region. 

It is easy to see that the two-particle current persists also above the gap voltage, eV > 2A: when the node 
«o is inside the gap, \E\ < A, the subgap circuit should consist of two resistors no matter how large the applied 
voltage is. Furthermore, since the singular point E = A always belongs to the integration region, the current will 
formally diverge at all voltages eV > 2A; generally, the n-particle current (n > 1) taken in the SPT approximation 
diverges at all voltages above eV = 2 A/ (n — 1). This catastrophe is known since the pioneering calculations of the 
two-particle current within the tunnelling Hamiltonian model 1 1 1 but it can be eliminated by using the improved 
perturbation expansion of section 15.21 This implies in fact that the currents have nonanalytical dependencies on 
the tunnelling parameter W at large voltages. 

We also note that the three-particle current disappears at large enough voltage, in contrast to the two-particle 
current which persists at all voltages above A/e. This is obvious from the circuit geometry in figure [3] as soon as 
eV exceeds 2A, the network period becomes larger than the energy gap, therefore the subgap circuit may involve 
no more than two resistors. This is relevant for all n-particle currents with n > 2 which persist only within the 
voltage intervals 2A/n < eV < 2A/(n — 2) and abruptly disappear at larger voltages (see figure[4]i. 
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7.1. Two-particle current at eV > 2A 

The two-particle current at eV > 2A is given by equation 

I 2 = ±r^ = ±[ A d E^L. (56) 
eRJo po + pi e/?7o #i 

To evaluate the integral, we use the IPT equations (f39b— (f42b to calculate N, while the functions N± \ are taken in 
the BCS form. Furthermore all the smooth functions can be approximated with their values at the singular point 

E = A, 

g(eV) = ^ k=±l N s (A + keV), {51a) 

JVT 1 +ATi =Y tk=±1 Ns 1 (A+keV) = h- 1 (eV). (57b) 

As it was expected, the two-particle current in this region is given by a nonanalytical expression with respect to W, 
and it is larger than the current value at smaller voltages eV < A, 

h = ?-C 1 p{eV)h{eV)~W b , (58) 
eR 

9\/3/4, ID junction, 

V2j planar junction. 

Here p(eV) is given by d39l or d4Tb depending on the junction geometry. 

Important property of the tunnelling IVC is the excess current / eX c, i.e., voltage-independent deviation of the 
total current from the ohmic IVC at large voltage eV ^> A. The excess current is readily evaluated by considering 
this limit in (fTTTb and d58i 

ID junction, 

, . . (59) 
planar junction. 



Ci = / de Imz(e) = 




Our analysis is not complete yet because the current in ( I58t grows infinitely when the voltage approaches 
2A. This divergence, caused by the singularity of N-i in (I57al) . can be eliminated by using a more accurate 
approximation d45ll-(l48l for N. Neglecting a nonsingular term N\ in the denominator in d56i l and approximating 
N\ with iV(3A) = 3 /2y/2, we obtain 

, , 4 f A 6qA f°° , x 

7 2 (2A)w— / dENNi = -^- / deRez(-e) 




V2eRJo 
ID junction, 

, . . (60) 
planar junction. 

Thus we see that the two-particle current possesses a pronounced peak at eV = 2A, which exceeds not only the 
current threshold value, but also the large excess current, see figureH] 



7.2. Three-particle current at eV > A 

The three-particle current at eV > A is given by dBTl ). in which the integration interval is now eV — A < E < A. 
Using the symmetry relation N(E) = N(—E), we reduce the integration region to the interval eV/2 < E < A, 
containing only one dangerous point E = A, 

6 f A dE 

h = — ■ (61) 

eR JeV/2 p-i 4-po + Pi 

At this point, all terms in the denominator turn to zero being calculated within the SPT approximation. To eliminate 
the divergence, we apply the IPT approach, d42b and d39l ), to calculate N, which then acquires a finite value ~ W~ h 
at the singular point. In the present case the function g defined in d38l has a complex-valued form g(eV) = |g|e 1<p 
because the energy £_i in this equation occurs inside the gap, 



A 3 / 2 eV-A f2A + eV 

g = — — , tanc> = \ . (62) 

11 ^eV[(2A) 2 -(eV) 2 } eV + A\2A-eV 

Thus the scaling factor p in d39l , ( HH and d42l is to be defined through the modulus of g, while the phase factor 
will remain in the equation for z in the ID geometry, 

z 3 exp(2i<j!») + (zVe-l) 2 = 0, (63) 
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and in the expression for z in planar geometry, 

z= [e-iexp(i<p)]- 1/2 . (64) 

Therefore the normalized spectral function becomes dependent on eV. The function N-\ can then be evaluated 
from the SPT approximation, keeping only large quantity N in the rhs, 

N-i =W\A/£-i\ a N. (65) 

At the singular point, N-i ~ W l ~ b \ the other relevant functions N\ and N-2 can be taken in the BCS form. 
Retaining only the largest resistance p_ i ~ W b ~ 1 in doTl l, we get 



6 f A 6WpA 
h = — dEN- 1N-2 = — |-tf-2 
eR JeV/2 eR 



C 2 ~W 



\+b 



1-1 

C 2 (eV) = / d£lmz(e,eV) (66) 



(for the planar model, C2(eV) = 2cos(<p/2 + 7l/4)). Here the functions N-2 and |_i should be taken at the point 
E = A. This expression diverges at eV = A and eV = 2A. In fact, the current has a finite peak value at eV = A; 
another peak appears slightly below eV = 2A, because the current turns to zero at eV — 2A, as follows from doTI ) 
(see figure|4|. 

At eV = A, the function N-2 is also large at the point E = A and should be kept in d65l ) together with N, 
i.e., N-i = W(N + N-2)', both these functions are to be evaluated using the IPT scheme d38~l>- (l4"2"l >. This leads to 
expressions Af = (2p)~ l Imz+(e) andAf_ 2 = Rez_(e), where the functions z± are given by the solutions of 

algebraic equations 

z 3 + + (z+Ve-l) 2 = 0, i£. + (z_Ve-l) 2 = (67) 

in the case of ID junction, or by explicit expressions 

z± = (£Tir 1/2 (68) 

for a planar junction; the scaling parameter p is (W 2 /6) 1//3 and (W/v / 3) 1/ ' 2 , respectively. Then the largest 
resistances are po ~ P-i ~ W , and we obtain 

6W r A 3WA f°° 

7 3 (A) = — / dENN-2 = / de Imz+Rez_ 

eR J a/2 eR Jo 

WA f 3.16, ID junction; 

= x < (69) 

eR I 37Z/4, planar junction. 

Analysis of the current behaviour near eV — 2A, which we do not present here, gives the following estimate 
for the current peak value: h max ~ W 4 / 5 A/eR for the ID junction and Ij, max ~ WA/eR for the planar junction. 



7.3. Four-particle current at eV > 2A/3 

The four-particle current at eV > 2A/3 is given by d54l ), where the upper integration limit is replaced by the value 
A — eV, representing a singular point of the integrand. The functions Af and A^_i are calculated from the SPT 
equations ( f3Tb , neglecting their values in the rhs, 

N = W\A/^\ a N u N-i =W|A/|_i| fl AL 2 . (70) 

Within the voltage interval 2A/3 < eV < A we can apply the BCS approximation for the functions N±2. 
However, the function N\ must be calculated by means of the IPT d38l- (l42l . as soon as the energy E\ hits the 
singular point. The resulting equations for N\ coincide with (I62l-d64li for in the previous section. Thus, at the 
singular point, A^i ~ W~ h , No ~ W l ~ h , N-\ ~ W, and N±2 ~ 1. Therefore the resistance po ~ W h ~ 2 dominates, 
which gives 



8AW 2 

h = —^-C 2 (eV)p(eV) 
eR 



A 2 



Ar_ 2 ^W 2+fe , (71) 



where the functions and N-2 are taken at E = A — eV. This expression diverges at the points eV —2A/3 

and eV = A, where the IPT should be applied. Since the SPT approximation for N-2 diverges at eV = 2A/3, we 
evaluate both the functions N\ and N-2 in the rhs of ( T70b using the IPT equations d38l-(l42l. 

N l (E) = (2p)- l lmz+(e), N- 2 {E) = (2p)- 1 Re Z -(e). (72) 
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The functions z± obey the equations d67| i or d68l l modified by the phase factors similar to d63l and d64l ). In this 
case, tan<p = — V2/5, and the scaling factor is p = |(W 2 /2) 1,/3 for ID junctions and p = (3WV3/8) 1 / 2 forplanar 
junctions. The corresponding estimates for the DOS functions are N\ ~ A^_2 ~ W~ b , N ~ iV_i ~ W l ~ b , while 
A^2 ~ 1 can still be taken in the BCS approximation. Therefore the largest resistance is po ~ W 2b ~ 2 , and we get 

4AW 2 /9V r 
*(2A/3)= — (-j ^ d £ Imz + Rez_ 

AW 2 [6.63, ID junction; 
= — ;r~ x < (73) 
eR I 5.26, planar junction. 

Analysis of the current peak near eV = A gives the estimate hmax ~ ffA/efl. 



8. Effect of first harmonics 



In this Section we evaluate the effect of higher harmonics on the dc subgap current. We restrict our calculation 
to first order in the tunnelling parameter W thus taking into account only two first harmonics m = ±1. The main 
conclusion of our calculation will be that including harmonics produces just insignificant quantitative changes, 
while major qualitative properties of the SGS, positions and scaling of the current steps will not change. We start 
with a general equation (TBI) for the dc current, and evaluate an additional contribution due to the first harmonics 
of the boundary values of the Keldysh and Green's functions, 
/•<*> fig 

51 = /— 32^ Tr ^=±l mt7z [M£>0)G K (E,m) 

+h(E,-m)G K (E,0)] (74) 
f°° AE 

= i / 4^L=±i m l V y( E > m ) Imv + V Imv,.(£,m)] . 

Here and below we use the subscripts (x,y,z) to indicate the matrix components of the first harmonics of the 
Green's and Keldysh functions, while the zero harmonics will be used as before without such subscripts. Thus v 
[V] in (l74l indicates y-component of g(E,Q) [G K (E,0)], and v y (E,±l) [V y (E ,±l)] indicates the y-component of 
£(E,±1)[G K (£,±1)]. 



8.1. Perturbation theory for the Green's functions 

To evaluate the current in d74l >. one needs to find the y-component in the Green's function expansion over the 
Pauli matrices, g(z,E,m) = G-u z + ia y v y + O x v x + w. In what follows, we will focus on the case of ID junction. 
Considering the Usadel equation ^ to first order in W, we arrive at the equation for the functions v y (z,E,m) and 
u z (z,E,m), 

meVu z = 2iA (u m u'' — v m v''j . (75) 

In this and following derivations we neglect small deviation of the functions u and v from the BCS form ( II 1M . We 
use the following convention: the zero harmonics with shifted energy arguments E + meV will be denoted with the 
subscript m, as before, e.g., u m = u(E + meV,0), while the first harmonics will be denoted with an argument, e.g., 
v y (E,±l); the absence of the explicit argument would mean the relevance to the both harmonics m = ±1. 

The function u z can be excluded from d75l l by virtue of the normalization condition («i + u~\)u z = (vi + 
v_i)v v following from (0. The boundary condition for d75l ) results from d9ol > and determines the boundary value 
of first derivative, 

v' y \z=-o = (W/2)v(l +mim_i +viv_i). (76) 
Solution of d75l l. v y (z) = v y exp(feiz) at z < 0, leads to the following relation at the boundary, 

v y (E ) m)=k^ V ' y (E ) m) ) fc 2 = + |_i)/2iA. (77) 

Equations (l76l l and (TTTl i allow us to express the first harmonics v y (E,m) through known zero harmonics u and v, 
and to establish the relation V y (E, 1) = v y (E, — 1) which implies that the second term in d74T > turns to zero. 
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8.2. Perturbation theory for Keldysh function 

To calculate the harmonics of the Keldysh function, it is convenient to separate the contributions of the harmonics 
of the Green's function and the distribution function, 

G K (£,w) = <5G K (£,m) + S K (£,m), 

8G K (E,m) =g R (E,m)ff m ~f m g A (E,m), (78) 

9 K (E,m)=glf(E,m)-f(E,m)g A m - 

In d78l ) we used the rule ([T2l > for the convolutions valid for the first harmonics, (AoB)(E,m) = A(E ,m)B- m + 
A m B(E,m), which allows us to write the equation for the function S K following from ( f8at and ( !%! > in a symbolic 
form 

[H, oS K ] = 2iA«9,(| R o d x § K + S K o d x g A ) , (79) 

where a small gradient of the distribution function, d x f ~ L 1 , has been neglected. The boundary condition to d79l ) 
reads 

| R o^§ K = W(| R oJ-Jo| A ), (80) 

where !T = g R o (/ — /) — (/ — /) o g A . According to ( T74b and (F78T >. we are interested in the y-component in 
the expansion § K = C7 Z U Z + io y V y + O x V x + W. Equation for this component, meVli z = 2iA (u m U!' — v m V'), 
follows from d79l l and looks similar to ( 1751 ). The normalization condition dHfcl i gives the relation (v m — v*_ m )'\l z = 
(u m — u*_ m )V y which allows us to obtain a closed differential equation for V y . Solution V y (z) = V y exp(q\z) of this 
equation at z < gives the following relation at the boundary, 

V y (E,m) = q[ l V y {E,m), q\ = & - <^*i)/2iA, (81) 

where the quantity V' is to be determined from the boundary condition ( |80l , 

V;| z =-o = W (1 - - v m v*_ m ) 4>(£,m), (82) 



<K(£,m)=M(/-i £ / fc ) + iM s sgn(m) £ 

i=±l i=±l 

where M s = Imv. At T = we obtain 4>(£,m) = 0(eV- |J?|)[Msgn(J?) +iMs sgn(m)]. 
8.3. Current of first harmonics 

The function V y (E,m), which determines the current in (1741 , is now expressed, according to ( F78b . through the 
quantities calculated in d76l)-d77l) and (|8T1l-([82l, V y (E,m) — v v /_,„ + v* f m + V y . Collecting all these equations and 
substituting them into d74b . we get 



51 = Zr II dElmv i - /- 1 ) Im ( vk i 1 cosh2 ^ ) 

+ Im [v^ 1 (/ - f-i) + v*q\ 1 (/ - h ) sinh 2 x] } , 
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Figure 6. 1-V characteristic with tunnelling SGS for planar junction computed numerically for the tunnelling parameter 
W = 10~ 3 . The inset shows a spike of the IVC near the threshold eV = 2A/3 of the three-particle current. 



Table 1. Threshold and peak values of the normalized multiparticle currents I„eR/A, n = 1 -4-4. The left and right 
sub-columns correspond to the ID and planar models of the tunnel junction, respectively. 

n eV = 2A/n eV = 2A/(n-\) eV = 2A/(n-2) 

Ik/2 — — 

2 nW 2.32W 2 / 5 2.50W 1 / 3 6.19W 1/3 2.83W 1/2 

3 6.33W 2 6.71W 2 3.16W 2.36W 0.44W 4 / 5 3.57W 

4 12.9W 3 14.9W 3 6.63W 2 5.26W 2 0.57W 0.48W 



where X = \{9l + &-\) an d X = + At T = 0, this equation simplifies: 

2W f eV 

8l=—j^ dE Im vim [v^f 1 cosh 2 ^ + ^ 1 sinh 2 ^)] . (83) 

Similar considerations in the case of a planar junction result in replacements q\ — > q\ and k\ — > k\ in the 
rhs of d83l . Noting that at eV < A the energy £_i in ( f83b appears in the subgap region, where 91 l = 0-i +in 
and §*j = — we find that 5/ turns to zero at eV < A, similar to /2. Numerical calculations show that the 
contribution of the first harmonics to the net dc current does not exceed 30% (see figure|5). From this we conclude 
that the adopted quasi-static approximation, where the nonzero harmonics are neglected, gives a relatively good 
approximation to a complete solution. 



9. Discussion 



Our analysis of the high-order multiparticle currents shows that they exhibit a similar pattern of the voltage 
dependence (see figure |4j: an n-particle current appears above the threshold eV = 2A/n, having roughly the 
value /„ - {2W) n ~ l h, then increases and shows dramatic peak while approaching eV — 2A/(n— 1); then it 
remains anomalously large within the voltage interval 2A/(n — 1) < eV < 2A/ (« — 2) and eventually disappears at 
eV = 2 A/ (n — 2), showing another strong peak at slightly smaller voltage. For convenience, all the threshold and 
peak values of the first four currents are brought together in the table Q] 

The net tunnel current consists of the sum of the n-particle currents, and therefore exhibits a pronounced 
step-like structure on the IVC with steps at the gap subharmonics eV = 2A/n, as shown in figure [6] obtained by 
numerical calculation at T = 0. The peaks of the multiparticle currents with numbers n + 1 and n + 2 produce 
small spikes at the nth threshold with n > 1; the example of such spike is presented in the inset of figure [6] The 
numerical procedure involves solving the set of recurrences d34t or ( f35l > for the functions 0^ which correspond to 
the subgap energies, < A (— < k < N+); the nonsingular terms in these equations are replaced with the 
BCS values. For the voltage values equal to the gap subharmonics one more equation is to be added as explained 
in the previous sections. We note that the results for the ID and planar geometries differ insignificantly for equal 
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values of the tunnelling parameter ( fTOb and (|28l i; in logarithmic scale, the difference can be detected only in the 
immediate vicinity of the peaks. 

Such a picture is quite similar to the tunnelling SGS in quantum point contacts J2] [3] |4] Qj] [T3], and the 
resulting IVC is found to be very close to the result for a point contact with the effective transparency D e s = 4W '. 
Thus, according to the definitions ( fTOb and (f28j) of the tunnelling parameter W, the enhancement factor D e g/D for 
the SGS scaling is equal to 3^q/£ in the ID junction and 3^/di in the planar junction. In particular, for planar Al 
junctions with £ ~ d = 50 nm and <^o = 300 nm, the enhancement factor may approach the value 100. 

The fact that the SGS in planar junctions is sensitive to the properties of the junction electrodes has important 
implications for characterization of the junction tunnelling layer. Indeed, the thickness of this layer in realistic 
junctions is inhomogeneous: there are spots with enhanced transparency, which mostly contribute to the tunnel 
current. If the linear sizes of such spots are large compared to the electron mean free path in the electrodes (in 
practice, the thickness of thin-film electrodes), the junction can be considered as a quasi-planar one, and the SGS 
should be enhanced according to our theory and depend on the electrode thickness. However, if such spots are 
small compared to the electron mean free path, the tunnel current rapidly spreads out in immediate vicinity of the 
spot without being affected by the impurity scattering, as it is in the ballistic constrictions; in this case, there must 
be no dependence of the SGS on the electrode properties in accord with the mesoscopic theory prediction dOH]. 
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